#include <stdio.h>
#include <sys/time.h>

#define POS(i, j, width) i * width + j
#define uc unsigned char

typedef struct {
	uc r, g, b;
} PIXEL;

// function to check errors CUDA functions
void CUDA_CHECK( cudaError_t error ) {
	if( error != cudaSuccess ) {
		fprintf(stderr,"ERROR: %s\n", cudaGetErrorString(error) );
		exit(1);
	}
}

__global__ 
void smooth(PIXEL * in, PIXEL * out, int width, int height, int SMOOTH_SQUARE)
{
	int i = blockIdx.y * blockDim.y + threadIdx.y;
	int j = blockIdx.x * blockDim.x + threadIdx.x;

	if (i >= height || j >= width) return;

	// sums for each component of RGB 
	int sumr, sumg, sumb;
	sumr = sumg = sumb = 0;

	// considers 5x5 square 
	int k, l;
	for(k = i - SMOOTH_SQUARE/2; k <= i + SMOOTH_SQUARE/2; k++) {
		for(l = j - SMOOTH_SQUARE/2; l <= j + SMOOTH_SQUARE/2; l++) {

			// assumes zero for border pixels
			if( k < 0 || k >= height || l < 0 || l >= width )
				continue;

			// acumulate
			sumr += in[POS(k,l,width)].r;
			sumg += in[POS(k,l,width)].g;
			sumb += in[POS(k,l,width)].b;
		}
	}

	// calculate the mean
	out[POS(i,j,width)].r = sumr/(SMOOTH_SQUARE * SMOOTH_SQUARE);
	out[POS(i,j,width)].g = sumg/(SMOOTH_SQUARE * SMOOTH_SQUARE);
	out[POS(i,j,width)].b = sumb/(SMOOTH_SQUARE * SMOOTH_SQUARE);
}

int main(int argc, char ** argv)
{
	struct timeval start_time_total, end_time_total;
	gettimeofday(&start_time_total, NULL);

	// require arguments to run
	if( argc != 4 ) {
		printf("usage: ./par-cuda <file_input_name.ppm> <file_output_name.ppm> <size_of_square_side(odd-number)>\n");
		exit(1);
	}

	// open file
	FILE * fp = fopen(argv[1], "r");
	FILE * fout = fopen(argv[2], "w");
	char encoding[4];
	uc maxValue;
	int width, height;
	int i, j;
	int SMOOTH_SQUARE = atoi(argv[3]);

	// read file header
	fscanf(fp, "%s", encoding);
	fscanf(fp, "%d %d", &width, &height);
	fscanf(fp, "%hhu", &maxValue);

	// allocate memory in host	
	PIXEL * in = (PIXEL *) malloc ( width * height * sizeof ( PIXEL ) );

	// read image
	for(i = 0; i < width * height ; i++) {
		fscanf(fp, "%hhu %hhu %hhu", &in[i].r, &in[i].g, &in[i].b);
	}

	struct timeval start_time_mem, end_time_mem;
	gettimeofday(&start_time_mem, NULL);

	// allocate memory in device
	PIXEL * dev_in;
	PIXEL * dev_out;
	CUDA_CHECK(cudaMalloc( &dev_in, width * height * sizeof ( PIXEL ) ));
	CUDA_CHECK(cudaMalloc( &dev_out, width * height * sizeof ( PIXEL ) ));

	// copy the image from host to device
	CUDA_CHECK(cudaMemcpy( dev_in, in, width * height * sizeof (PIXEL), cudaMemcpyHostToDevice ) );

	// define the thread hierarchy
	dim3 threadsPerBlock(32, 32);
	dim3 numBlocks((width + threadsPerBlock.x - 1) / threadsPerBlock.x, 
		(height + threadsPerBlock.y - 1 )/ threadsPerBlock.y);

	struct timeval start_time, end_time;
	gettimeofday(&start_time, NULL);

	// call the kernel and wait for all thread finalize
	smooth <<< numBlocks, threadsPerBlock >>> (dev_in, dev_out, width, height, SMOOTH_SQUARE);
	CUDA_CHECK(cudaThreadSynchronize());

	gettimeofday(&end_time, NULL);

	// copy the result image to host 
	CUDA_CHECK(cudaMemcpy(in, dev_out, width * height * sizeof ( PIXEL ), cudaMemcpyDeviceToHost ));

	gettimeofday(&end_time_mem, NULL);

	// write the results in output file
	fprintf(fout, "%s\n", encoding);
	fprintf(fout, "%d %d\n", width, height);
	fprintf(fout, "%hhu\n", maxValue);
	for(i = 0; i < height; i++) {
		for(j = 0; j < width; j++) 
			fprintf(fout, "%hhu %hhu %hhu ", in[POS(i,j,width)].r, in[POS(i,j,width)].g, in[POS(i,j,width)].b);
		fprintf(fout, "\n");
	}

	free(in);

	cudaFree(dev_in);
	cudaFree(dev_out);
	fclose(fp);
	fclose(fout);

	gettimeofday(&end_time_total, NULL);

	//double result = end_time.tv_sec - start_time.tv_sec +  (end_time.tv_usec - start_time.tv_usec)/1000000.0;
	double result_mem = end_time_mem.tv_sec - start_time_mem.tv_sec +  (end_time_mem.tv_usec - start_time_mem.tv_usec)/1000000.0;
	//double result_total = end_time_total.tv_sec - start_time_total.tv_sec +  (end_time_total.tv_usec - start_time_total.tv_usec)/1000000.0;

	//printf("Only GPU processing: %lfs\n", result );
	printf("%lf\n", result_mem );
	//printf("GPU processing + memory copy + in/out - TOTAL : %lfs\n", result_total);

	return 0;
}
